Linear Stability Analysis of Differentially Rotating Polytropes 
— New results for the m = 2 /-mode dynamical instability — 
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ABSTRACT 

We have studied the /-mode oscillations of differentially rotating polytropes by making use 
of the linear stability analysis. We found that the critical values of T/|W| where the dynamical 
instability against the m = 2 /-mode oscillations sets in decrease down to T/|W| ~ 0.20 as the 
degree of differential rotation becomes higher. Here to is an azimuthal mode number and T and 
W are the rotational energy and the gravitational potential energy, respectively. This tendency 
is almost independent of the compressibility of the polytropes. These are the first exact results 
of the linear stability analysis for the occurrence of the dynamical instability against the m = 2 
/-modes. 

Subject headings: gravitational waves — stars: neutron — stars: rotation — stars: oscillation 



1. Introduction 

It is well known that rapidly rotating stars suf- 
fer from dynamical instability against nonaxisym- 
metric oscillations. In particular, the growth rate 
for the to = 2 or bar-type /-mode instability 
has been considered to be largest. The classi- 
cal results about this instability were first stud- 
ied by Riemann (1860). Riemann found that uni- 
formly rotating incompressible fluids become dy- 
namically unstable when the ratio of the rota- 
tional energy, T, to the absolute value of the grav- 
itational potential energy, W, exceeds 0.27, i.e. 
[3 = T/\W\ > 0.27 (see also Bryan 1889). 

Around 1970, Ostriker and coworkers developed 
a new numerical method to solve equilibrium con- 
figurations of rapidly rotating compressible stars 
(Ostriker & Mark 1968; Ostriker & Bodenheimer 
1968; Mark 1968; Jackson 1970; Bodenheimer & 
Ostriker 1970, 1973; Bodenheimer 1971). Ostriker 
and coworkers also developed a new technique to 
explore stabilities of rotating equilibrium configu- 
rations by extending the tensor-virial technique of 
Chandrasekhar ( Chandrasekhar 1969; Tassoul & 
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Ostriker 1968; Ostriker & Tassoul 1969; Ostriker 
& Bodenheimer 1973). Their main result about 
the to = 2 dynamical instabilities of compressible 
stars was that rapidly rotating configurations be- 
come unstable when the condition (3 > 0.27 is sat- 
isfied irrespective of the compressibilities and the 
rotation laws as far as they investigated. Although 
there was no proof of this "universal" nature of 
the critical value /3 C ritical ~ 0.27, most people have 
believed that the critical value for the m = 2 /- 
mode dynamical instability is 0.27. This "belief 
has lasted rather long, although there appeared a 
paper that shows the tensor virial technique for 
stability analysis for differentially rotating stars 
gives neither the necessary condition nor the suf- 
ficient condition for the stability (Bardeen et al. 
1977). 

However, the situation has been changing be- 
cause there appear many recent results of hydro- 
dynamical simulations which show that configura- 
tions with lower values of fi suffer from dynami- 
cal instabilities (see e.g. Tohline & Hachisu 1990; 
Pickett et al. 1996; Toman et al. 1998; Brown 
2000; Centrella et al. 2001; for simulations of lin- 
earized equations see e.g. Liu 2002). Rough values 
of the critical configurations against nonaxisym- 
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metric perturbations are /3 C riticai ~ 0.16 (to = 2) 
for self-gravitating toroidal polytropes with the 
polytropic index N = 1.5 (Tohline & Hachisu 
1990), /3 cri ticai ~ 0.2 (to = 2) for N = 1.5 
polytropic star-disk like configurations (Pickett et 
al. 1996), Critical ~ 0.24 - 0.28 (to = 2) for 
N = 5/4,6/4,10/4 polytropes with the rotation 
laws that the distribution of the specific angu- 
lar momentum is that of uniformly rotating poly- 
tropes with the same polytropic indices (Toman 
ct al.1998), Critical ~ 0.24 (m = 2) for N = 1.5 
polytropes with the angular velocity of the expo- 
nential decay with the distance from the rotation 
axis (Brown 2000), and /^critical ~ 0.14 (to = 1) 
for toroidal-like N = 3.33 polytropes (Centrella et 
al. 2001). Liu found /3 crit icai ~ 0.25 (to = 2) for 
the angular momentum distribution which could 
be realized after accretion induced collapse of uni- 
formly rotating massive white dwarfs, although 
his simulations were carried out by using the lin- 
earized equations. 

These results seem to imply that the critical 
value of (3 for the occurrence of the to = 2 bar- type 
dynamical instability cannot be considered to be 
universal any more. In order to see the unstable 
nature of rotating stars, we need to investigate se- 
quences of rotating configurations systematically. 
For that purpose, the best method would be the 
linear stability analysis rather than hydrodynam- 
ical simulations. 

Except for non-self-gravitating disks, linear sta- 
bility analyses for rapidly rotating self-gravitating 
stars have been done only for limited kinds of com- 
pressibilities and rotation laws, and for limited 
types of oscillations. Luyten (1990, 1991) found 
a new dynamical instability for differentially ro- 
tating polytropes. This instability seems to be- 
long to the shear instability which is also found 
in the theory of accretion disks (see e.g. Pa- 
paloizou & Pringle 1984; Narayan & Goodman 
1989). This instability appears for differentially 
and very slowly rotating stars and disappears for 
differentially but very rapidly rotating stars. How- 
ever, Luyten's analysis was applied only to very 
limited kinds of rotation laws and the displace- 
ments were restricted only in the plane parallel to 
the equator. Very recently, Shibata et al. (2002, 
2003) have investigated equilibrium sequences of 
differentially rotating polytropes by using linear 
stability analysis, the numerical scheme of which 



is the same as that employed in this paper, and 
shown that the dynamical instability of a new type 
can set in even when (3 is ~ 0.05 if the degree of 
differential rotation is very high. 

However, even the classical m — 2 /-mode in- 
stability has not been fully studied by the linear 
stability analysis. Thus, in this paper, we will 
study systematically the dynamical stability of 
rapidly and differentially rotating polytropic stars 
by using the linear stability analysis. It should 
be noted that the linear stability analysis of dif- 
ferentially rotating polytropes in this paper will 
give exact values of /3 C riticai for the onset of m = 2 
/-mode dynamical instabilities for the first time. 

2. Numerical Method 

The method to investigate stability by solving 
linearized basic equations is the same as that used 
in Karino ct al. (2000, 2001). Therefore, we will 
summarize the essence of the method briefly. 

2.1. Equilibrium states 

Axisymmetric equilibrium configurations of ro- 
tating Newtonian stars are obtained by solving the 
following equations in the spherical polar coordi- 
nates (r, 6, <p): 

^-Vp= -\7(f> + r sin 9n 2 e R 7 (1) 

Acf) = AirGp , (2) 

where p, p, <fi, ft, e*R, and G are the pressure, the 
density, the gravitational potential, the angular 
velocity, a unit vector in the R direction of the 
cylindrical coordinates (R,ip,z), and the gravita- 
tional constant, respectively. As for the equation 
of state, we assume the polytropic relation as fol- 
lows: 

P = Kp 1+ ^ , (3) 

where K and N are the polytropic constant and 
the polytropic index, respectively. The rotation 
law is assumed as follows: 

n c A 2 

° = (r sin 6Y+ A 2 ' (4) 

where Q c is the central angular velocity. The 
quantity A is a parameter which represents the 
degree of differential rotation. The rotation be- 
comes more differential as A becomes smaller. We 



2 



call this rotation law as the j-constant rotation 
law. 

One equilibrium state can be characterized by 
specifying three parameters as follows: the poly- 
tropic index N, the parameter A and the param- 
eter describing the amount of rotation which can 
be represented by the following axis ratio: 

q = r p /r e , (5) 

where r p and r c are the polar and equatorial radii, 
respectively. These equations are solved by ap- 
plying the proper boundary conditions and using 
the SFNR method of Eriguchi and Muller (1985). 
Equilibrium models in this paper have been solved 
by using mesh numbers (r x 6) = (42 x 19). 

2.2. Perturbed States 

As mentioned before, the linearized equations 
are solved by the scheme developed by Karino et 
al. (2000, 2001). The perturbed quantities are 
expanded as follows: 

Sf{r,6,ip,t) =^2exp(-i(at-mtp))f m (r,6), (6) 

m 

where Sf means the Euler perturbation of the cor- 
responding quantity. Here, m is the azimuthal 
mode number and a is the eigenfrequency of the 
oscillation mode. As for the density and the 
pressure perturbations, we assume the adiabatic 
changes. Since we want to study the bar-type 
oscillations, we will consider only the oscillation 
modes with m — 2 in this paper. 

In the stability analysis in this paper we have 
also used mesh numbers (r x 9) = (42 x 19). Al- 
though this mesh number seems to be too small, 
the obtained eigenvalues can be considered to be 
accurate enough. We have checked the eigenvalues 
by changing the mesh number. For N — 1.0 poly- 
tropes with extremely rapid rotation (q = 0.2), the 
accuracy of the obtained eigenfrequency can be 
estimated by 5a(m,n) = \(a m — <J n )/a m \, where 
(Tfe denotes the eigenvalue when the mesh num- 
ber in the r-direction is k. For A = 1.0 models, 
6er(42,30) = 1.8 %, and <5er(42,36) = 0.6 %. 

3. Numerical Results 

By fixing the values of A and N and by chang- 
ing the value of q, we can obtain eigenvalues of 



m = 2 /-mode oscillations along an equilibrium 
sequence. After that by changing the values of 
A and/or N and following the same procedure, 
we can obtain eigenstates for many equilibrium 
sequences. In this paper, we have examined se- 
quences with N = 0.0,1.0 and 1.5, and mainly 
1.0 < A^ 1 . It may be noted that for models with 
large values of A, mass begins to shed from the 
equatorial surface before the dynamical instabil- 
ity sets in. Hence it is not necessary to consider 
larger values of A more than the value for which 
mass shedding occurs. 

In Figure 1, the real part of eigenfrcquencies 
of the N = 1.0 sequences with A = 1.0 is plot- 
ted against the value of 0. In this plane, the 
m = 2 /-mode has two branches for lower values 
of 0. The upper branch corresponds to the co- 
rotating mode and the lower one corresponds to 
the counter-rotating mode. These two branches 
merge into two complex conjugate branches for 
higher values of 0. The positive imaginary part is 
shown in Figure 2. The model for which the two 
branches merge corresponds to the critical config- 
uration from which dynamically unstable config- 
urations follow. For this sequence, i.e. N = 1.0 
and A = 1.0 sequence, the critical value where the 
dynamical instability for m = 2 /-mode sets in is 
~ 0.266. 

The eigenvalues of configurations of N = 1 
polytropic sequences with several values of A are 
shown against the value of 0: real parts in Fig- 
ure 3 and imaginary parts in Figure 4. It is clear 
from Figures 3 and 4 that the values of where 
two real modes change to two complex conjugate 
modes become smaller as the degree of differen- 
tial rotation becomes higher. It means that the 
critical point of dynamical instability will appear 
for smaller values of as the degree of differential 
rotation becomes significant. 

The dependency of the critical points of dynam- 
ical instability on the degree of differential rotation 
and the compressibility can clearly be seen from 
Figure 5. In this figure, the values of at the 
critical points for N = 0.0, 1.0 and 1.5 polytropes 
are plotted against the parameter A^ 1 . The crit- 
ical values of for the m = 2 /-mode dynamical 
instability depend significantly on the values of A 
but depend only slightly on the values of N for the 
high degree of differential rotation. Note that, as it 
should be, the critical value of /^critical approaches 
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to /^critical = 0.27 when the rotation law tends to 
the rigid rotation (A^ 1 — ► 0) and the unperturbed 
configuration becomes the Maclaurin spheroid. 

In Figures 6-8, the eigenfunctions of the density 
perturbation for N = 1 polytropes with A = 1.0 
but with different values of q are plotted against 
the normalized distance from the center of the star 
at selected values of 0. In slowly rotating regions, 
the eigenfunctions are almost similar to those of 
/-mode oscillations for the Maclaurin spheroid, 
i.e. monotonically increasing functions of the dis- 
tance from the center (see e.g. Tassoul 1978). As 
the stellar rotation becomes faster, however, the 
eigenfunctions on and near the equator attain their 
maximum values on the way to the surface but not 
on the surface. 

4. Discussions 

By making use of the linear stability analysis, 
we have shown that, for m = 2 /-mode oscilla- 
tions, higher degree of differential rotation makes 
the stars more unstable than previously thought. 
Although results of recent dynamical simulations 
have shown the same tendency, it should be noted 
that our present results have been obtained by the 
systematic analysis of the oscillations and that we 
can understand the occurrence of this type of in- 
stability from a more rigorous study. 

Since the parameter regions for the existence of 
dynamically unstable configurations are widened, 
we need to reconsider some astrophysical situa- 
tions under the light of the present results. There 
are several situations where high differential ro- 
tation would be realized. One example can be 
outcomes of collapsing gases or stars, i.e. star 
formation stages or core collapse stages of mas- 
sive stars or accretion-induced collapse processes 
of stars. Other example may be outcomes of col- 
lision or merging of two stars, i.e. merging of two 
neutron stars in a binary neutron star system. 

In the star formation process, if the angular mo- 
mentum is not lost from the gaseous system, the 
rotation of the final young stellar objects would 
be extremely rapid and highly differential. Thus 
they might suffer from the dynamical instability 
discussed in this paper. Although the outcomes 
of the dynamical instability have not been fully 
understood yet, several authors have argued that 
dynamical instabilities during the star formation 



stages are related to binary formations (see e.g. 
Bonnell 1994; Matsumoto & Hanawa 1999), and 
others have suggested that they would result in 
new systems consisting of a central young stellar 
object and a disk around it (see e.g. Bate 1998). 
In order to obtain a final answer to this problem, 
reliable nonlinear simulations and high-resolution 
infrared/radio observations about star forming re- 
gions are required. 

Another example is related to newly-born neu- 
tron stars. There are several ways for neutron 
stars to be formed. One is the formation via core 
collapse of massive stars followed by Type II su- 
pernova explosions. In recent simulations, it is 
shown that the compact object formed from an 
iron core collapse is rotating rapidly and differ- 
entially with (3 ~ 0.2 (Zwerger & Miiller 1997; 
Dimmclmcier et al. 2002). Another is a merg- 
ing of two neutron stars in a binary neutron star 
system. In the merging process in a binary neu- 
tron star system, there appears a rapidly rotating 
massive neutron star with high degree of differ- 
ential rotation (see e.g. Shibata and Uryu 2000). 
Neutron stars are also considered to be formed by 
accretion onto massive white dwarfs. In this ac- 
cretion induced scenario, the outcome would be 
a rapidly and differentially rotating massive neu- 
tron star (see e.g. Liu & Lindblom 2000). There- 
fore, for these differentially rotating neutron stars, 
it would be very likely that the dynamical insta- 
bility would set in easily and lead the systems 
to nonaxisymmetric configurations just depend- 
ing on the degree of differential rotation and ir- 
respective of the compressibility. Those deformed 
and rapidly rotating compact objects would emit 
a large amount of gravitational waves. Therefore, 
such young neutron stars could be important tar- 
gets of ground-based laser interferometric detec- 
tors of gravitational waves (see e.g. Brown 2000; 
Shibata et al. 2002) 2 . 

Concerning the rotation law dependence of the 
critical values for the m = 2 /-mode dynamical in- 
stability, we have performed linear stability anal- 
ysis for equilibrium models with the v-constant 
rotation law (see e.g. Karino et al. 2001). How- 

2 Whcn wc consider these compact stars, the effect of general 
relativity becomes important. In relation to the dynamical 
instability, however, it is only argued that the critical values 
might become slightly smaller by the general relativistic 
effect (Shibata, Baumgartc, & Shapiro 2000). 
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ever, the effect of differential rotation seems very 
weak and the dynamical instability occurs only for 
models with very stiff equations of state and/or 
extremely small values of the parameter A. Ad- 
ditionally, the change of the critical values of /3 
is very little. The critical values for N = 0.5 
polytropes are T/\W\ ~ 0.266 for A = 0.8 and 
T/\W\ ~ 0.260 for A = 0.5. 

5. Summary 

In this paper, we have performed linear stability 
analysis of equilibrium sequences of differentially 
rotating polytropes in Newtonian gravity and ob- 
tained critical values of j3, i.e. /^critical, which de- 
termine the equilibrium states where the dynam- 
ical instability against m = 2 /-mode oscillations 
sets in. The rotation law which we have employed 
is the j-constant law (Eq. (4)). Our important 
finding is that the critical points of m = 2 /- 
mode dynamical instability depend strongly on the 
degree of differential rotation while they depend 
weakly on the compressibility. The critical values 
of /3criticai decreases as the degree of differential ro- 
tation becomes higher. Therefore, we should not 
use the "universal" value of /3 C riticai = 0.27 for the 
stability analysis of rapidly rotating objects with 
high degree of differential rotation from now on. 
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Fig. 1 . — Real part of the eigenvalue of the m = 2 
/-mode in units of \/AirGp c is plotted against the 
ratio T/|W| for polytropes with N = 1.0 and 
A = 1.0. Here p c denotes the central density. In 
slowly rotating regions, the /-mode has two real 
branches, but when the stellar rotation exceeds a 
certain critical point, these two branches merge 
into one branch, in this real part of the eigen- 
value plane, since eigenvalues arc complex conju- 
gate. This corresponds to the critical point where 
the m = 2 /-mode dynamical instability sets in. 
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Fig. 2. — Positive imaginary part of the eigen- 
value for the same mode as Figure 1 is plotted 
against T/|W|. In slowly rotating regions, the 
imaginary part of the eigenfrcquency vanishes and 
it implies that the configurations are dynamically 
stable. After the stellar rotation reaches a certain 
critical point, the imaginary part of the eigenval- 
ues appears and grows very rapidly. This corre- 
sponds to the critical point where the m = 2 f- 
modc dynamical instability sets in. 
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Fig. 3. — Real part of the eigenvalues around the 
merging points is plotted against T/|VF| for poly- 
tropic stars with N = 1.0 and A = 0.7, 0.8, 0.9, 1.0. 
As the degree of differential rotation becomes 
higher, the values of T/|W| at the critical points 
become lower. 
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Fig. 4. — Same as Figure 3 but for the imaginary 
part of eigenvalues. 
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Fig. 5. — Values of TYITF] at the critical points 
where the m = 2 /-mode dynamical instability 
sets in are plotted against the degree of differential 
rotation, A . The stellar models are polytropes 
with N = 0.0, 1.0 and 1.5. As the degree of differ- 
ential rotation becomes higher, the critical value 
of T"/|W| tends to decrease. 
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Fig. 6. — Eigenfunctions of the density perturba- 
tion for the bar-mode oscillation of a differentially 
rotating polytropic star with N = 1.0, A = 1.0, 
and q — 0.90 are plotted against the normalized 
distance from the rotation axis. Here R 8 (0) is the 
radius to the surface along the the polar angle 6. 
Three curves correspond to the eigenfunctions on 
the spokes with 9 — n/6 (curve 1), 6 = w/3 (curve 
2), and 9 = w/2 (curve 3). The value of the eigen- 
functions is normalized so that the maximum ab- 
solute value among all the perturbed quantities 
becomes unity. 
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7. — Same as Figure 6 but for q = 0.55. 



8. — Same as Figure 6 but for q = 0.30. 



